%% This file will generate bootstrap samples for:
% 1. The win probabilities of each firm
% 2. Win probabilities under counterfactual removing country-level fixed border costs
% 3. Win probabilities under counterfactual removing all border costs.

clear;

%setup file to keep track of iteration count...
logfile = fopen('boot_log.out','w');

%% Either run these scripts or load their output files
%setup; 
load data_structures;
%Main_mpec
load est_point;
%
% Now storing standard errors (and all stats) in est_point. 
%load est_se;

%name_vec = {'rhohat' , 'rho_nF' ,  'rho_nB' , 'natShare', 'natShareNF' , 'natShareNB' , 'markups' ,  'markupsNF' , ...
%    'markupsNB' , 'profits' ,  'profitsNF' , 'profitsNB' ,  'conSurp' , 'conSurpNF' , 'conSurpNB' , ...
%    'delasticity_ger' , 'delasticity_dnk' , 'fc_bounds' };


%%
%Draw B points from the distribution of beta0 
B = 200;
%K is pre-specified, but just in case
K = length(thetahat);

%Setup a random number stream
mystream = RandStream.create('mrg32k3a', 'seed', 110307);
RandStream.setDefaultStream(mystream);
draws = randn(K, B);

thetaB = bsxfun(@plus, thetahat, chol(Sig_hat, 'lower')*draws);



%% Setup Jac structure to solve for baseline rho 
%Setup Jacobian for model with fixed-costs...
%These are the sparsity templates for each project the first row, is the 
%fringe share and the final column is theadding up constraint
J_m_ger = [[ones(1,m.num_firms_ger);eye(m.num_firms_ger)],ones(m.num_firms_ger+1,1)];
J_m_dnk = [[ones(1,m.num_firms_dnk);eye(m.num_firms_dnk)],ones(m.num_firms_dnk+1,1)];

%This is the Jacobian of the Constraints by the market share parameters
%(rhos) (NOTE: the prime at the end of the definition!!!)
Jac_Pattern_FC = [kron(eye(m.num_pro_ger),J_m_ger), zeros((m.num_firms_ger+1)*m.num_pro_ger, (m.num_firms_dnk+1)*m.num_pro_dnk) ; %These are the german market constraints
    zeros((m.num_firms_dnk+1)*m.num_pro_dnk, (m.num_firms_ger+1)*m.num_pro_ger), kron(eye(m.num_pro_dnk), J_m_dnk) ]'; 


%r0 is the rho values that are in equilibrium for thetahat, NO LONGER USED
%

%m is the model specification with fixed costs.
%Set up bounds:
LB = zeros(size(rhohat));
UB = ones(size(rhohat));

ktropts = optimset('DerivativeCheck','off','Display','off',...
           'GradConstr','on','GradObj','on','TolCon',1E-6,'TolFun',1E-6,'TolX',1E-6,'JacobPattern',Jac_Pattern_FC);


%% The loop that solves the models and generates statistics


%if matlabpool('size') == 0
%   matlabpool(3);
%end

for var = 1:length(name_vec)
    eval(sprintf('%s_mat = [];', name_vec{var}));
end


for b = 1:B
    tic;
    % Solve model for baseline rho
    [rhohat, FVAL, EXITFLAG, OUTPUT] = ktrlink(@(x_0) dummy_objective(x_0), rhohat, [],[],[],[],LB,UB,@(x_0) model_constraints(x_0, thetaB(:,b),m),ktropts,'knitro.opt');
    if EXITFLAG
        disp('Warning: Exitflag baseline solve for rho is nonzero');
        fprintf(logfile, 'ITERATION %d Warning: Exitflag baseline solve for rho is nonzero\n', b);
    end
    
    
    % Do postestimation statistics
    
    [rho_nF rho_nB natShare natShareNF natShareNB markups markupsNF markupsNB ...
    profits profitsNF profitsNB conSurp conSurpNF conSurpNB ...
    delasticity_ger delasticity_dnk fc_bounds ] = postEstStats(m, m_nFC, thetaB(:,b), rhohat);
    
    for var = 1: length(name_vec)
        eval(sprintf('%s_mat = [%s_mat %s(:)];', name_vec{var}, name_vec{var}, name_vec{var}))
    end
    
    
    % Save stuff...
    save boot_chkpt;
    disp(sprintf('******ITERATION %d COMPLETE!!*******', b));
    fprintf(logfile, '******ITERATION %d COMPLETE!!*******\n', b);
    toc;
end

%ALERT: In the bootstrap we solve for "rhohat", while the estimate rho is
%saved under "e_rhohat" we copy it to the name "rhohat" to avoid conclusion
rhohat = e_rhohat;

%matlabpool('close');
fclose(logfile);
save counter_boot;


%% For each element of interest, compute standard error and reshape our variables into their standard format.
for var = 1: length(name_vec)
   eval(sprintf('%s_se = ((1/B)*sum((repmat(e_%s(:), 1, B) - %s_mat).^2,2)).^(1/2);', name_vec{var}, name_vec{var}, name_vec{var}));
   eval(sprintf('%s_se = reshape(%s_se, size(e_%s));', name_vec{var}, name_vec{var}, name_vec{var}));
   eval(sprintf('clear %s', name_vec{var}));
end


save counter_boot;
